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Summary 

The  present  iineil  Report  contadns  final  chapters  4,  5  and  6  describing  the  Computational 
Damage  Model  for  Thermoviscoelastic  Laminated  Composite  Materials.  The  detailed 
models  for  damageable  composite  elements  for  the  general  3D  case  was  described  in  [1], 
Chapters  1  and  2.  The  model  was  based  on  det^ed  physical  description  and  developed 
kinetic  equations  for  accumulations  of  damages  in  shear,  tension  and  delamination.  The 
model  constants  were  distinguished  that  shotdd  be  determined  from  the  experiments. 

The  Chapter  3  [2]  contained  the  description  of  the  experimental  procedure  on  twisting 
and  tension  of  thin-walled  tubidar  samples.  The  application  of  the  damage  model  for 
2D  cases  of  deforming  thin  shells  was  developed.  Some  solutions  were  obtained  within 
the  frames  of  quasistatical  approximation  enabling  to  develop  some  of  the  damage  con¬ 
stants  for  low  rates  of  loading.  Recommendations  for  developing  constants  by  comparing 
experimental  and  theoretical  results  were  worked  out. 

The  present  report  contains  the  derivation  of  the  model  equations  for  quasidynamical 
problems  on  twisting  and  tension  of  tubular  thin  shells  (Chapter  4).  The  solutions  of  the 
problem  are  obtained  for  different  rates  of  loading  (Chapter  5).  Based  on  the  obtained 
solutions  for  dynamical  behavior  of  parameters  under  different  rates  of  loading  (twisting 
and  tension)  detaiiled  recommendations  are  worked  out  for  the  test  procedures  wich  need  to 
be  performed  to  develop  the  constants  for  the  model  of  damageable  laminated  composite 
materi^ds. 

The  Chapter  6  contains  the  developed  mathematical  methods  to  model  numerically  the 
behavior  of  material  and  its  further  dynamical  deforming  after  accumulation  of  damages 
leads  to  destruction  of  material  in  the  damaged  zone. 


Chapter  4 


Quasidynamical  problems  of 
deforming  and  breakup  of 
thin-walled  cylindrical  tubular 
two-phase  laminated  composite 
samples 


We  reg£ird  the  classical  problems  of  quasidynamiccil  twisting  or  tension  of  a  thin-walled 
tubular  sample.  The  length  of  the  sample  is  denoted  i,  the  radius  (external)  a,  the 
thickness  of  the  walls  h  {h/a  <  1,  a/L  <  1).  Let  the  sample  be  loaded  monotonously 
by  the  growing  rotating  moment  M  =  M{t)  {M{t)  >  0)  applied  to  the  end  points  of  the 
sample  along  the  z-axis,  or  by  axial  tensile  forces  F  =  F{t)  {F{t)  >  0)  (Fig.  4.1). 

In  the  very  beginning  {t  =  to)  there  are  no  stresses  and  strains  in  the  sample  and  its 
temperature  is  uniform  (T  =  To).  Let  the  time  of  loading  until  breakup  be  much  larger 
than  the  characteristic  time  of  elastic  wave  traveling  along  the  sample  but  much  less  than 
the  characteristic  time  of  heat  propagation  along  the  sample  due  to  thermal  conductivity. 
Then  the  problem  could  be  regarded  within  the  frames  of  quasidynamical  approximation. 
In  fact  the  term  "quasidynamical”  means  that  we  apply  a  quasistatical  approximation  to 
mechanical  processes  and  adiabatic  approximation  to  the  thermal  processes. 


4.1  Twisting  of  a  thin- walled  tubular  sample 


In  case  of  pure  twisting  of  a  sample  F{t)  =  0  only  one  component  of  the  stresses  tensor 
(shear  stress  cgz)  differes  from  zero: 


<Tez{t)  = 


M{t) 

2'Ka?h' 


-4’ 


Let  (T0z{t)  =  f  •t',  f  =const>  0. 


h«a«l 


Then  the  system  of  constitutive  equations  for  the  problem  will  take  the  form: 
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L  1  —  a  zi 

The  constants  used  in  (4.1)  have  the  following  form  [2]: 

C22  -  C23  =  2  (/i) ,  a22  -  a23  =  — ;  .2  (-)  ? 
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K.=K-  (a  =  1, 2)  (4.2) 

where  Ea  -  isothermic  Young  modulus,  Aa  -  isothermic  Lame  coefficient.  The  physical 
meaning  of  all  the  constants  was  discussed  in  details  in  [1,2]. 


4.2  Tension  of  a  thin- walled  tubular  sample 


In  case  of  pure  tension  of  the  sample  {M{t)  =  0)  only  one  component  of  stress  tension 
(Tzz  differes  from  zero: 


(Tzzit)  = 


m 

2‘irah' 


Let  (Ttx{t)  =  f  •i,  f  =  const  >  0. 

Then  the  governing  system  of  equations  will  have  the  following  form: 
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The  constants  in  (4.3)  that  have  not  been  described  in  (4.2)  are  the  following: 
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Chapter  5 


Dyn2imical  behavior  of  governing 
parameters  in  deforming  and 
breakup  of  thin-walled  two-phase 
cylindrical  tubular  composite 
samples 


The  present  chapter  contains  the  results  of  numerical  investigations  of  the  damage  pa¬ 
rameters  behavior  for  the  model  laminated  composite  material  in  possible  experiments  on 
twisting  and  tension  of  tubular  samples  for  different  rates  of  loading.  Since  the  damage 
peirameters  cannot  be  measured  directly  in  the  experiment  the  data  on  the  behavior  of 
strains  and  temperatures  depending  on  the  damage  parameters  accumulation  is  also  en¬ 
closed.  The  results  are  mmed  to  illustrate  the  influence  of  occuring  damages  of  different 
type  (damages  in  sheeir,  tension  and  delamination)  on  the  behavior  of  sheeur  and  tensile 
strciins  under  different  rates  of  loading.  The  sharp  differences  in  strains  on  switching  on 
the  damage  parameters  make,  it  possible  to  determine  the  values  of  damage  pctrzuue- 
ters  [2]  ensuring  coincidence  of  theoretical  and  experimental  curves  for  shear  and  tensile 
strciins  under  different  rates  of  loading.  The  strain  history  is  given  for  each  rate  of  loading 
for  the  cases  accounting  for  and  neglecting  the  accumulation  of  damages,  illustrating  the 
differences  in  dynamicad  behavior  of  materials. 


5.1  Twisting  thin- walled  tubular  samples  at  differ¬ 
ent  rates 

In  pure  twisting  of  tubulcir  samples  by  a  rotating  moment  the  shear  stress  aez  was  assumed 
to  change  linearly  (Chapter  4): 

=  f 't. 
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(5.1) 


The  experimentally  measurable  values  e$z  and  T  under  these  conditions  can  depend  only 
on  djimages  in  shear 

(«) 

and  delamination  (wa)-  Now  we  will  illustrate  the  solution  of  the  direct  problem:  choosing 
the  model  composite  material  with  given  values  of  all  the  parameters  we  will  determine 
variations  of  strains  and  temperatures  under  different  rates  of  loading.  That  will  help  us  to 
determine  the  necessary  experiments  for  the  solution  of  the  inverse  problem  of  developing 
the  material  parcuneters. 

The  values  of  parameters  for  the  model  composite  material  are  given  in  the  table  5.1. 
AH  the  pareimeters  were  maintained  constant.  Only  the  twisting  stress  rate  /  (5.1)  varied 
from  10®  Pa/s  to  10^  Pa/s. 


Parameter 

Value 

Unit 

Description 

c 

0.61 

- 

matrix  (phase  1)  share  in  composite 

To 

300 

K 

initial  temperature 

pi 

1600 

kg/m® 

density  of  the  phase  1 

P2 

1300 

kg/m® 

density  of  the  phase  2 

Pi 

53  •  10^ 

Pa 

she2ir  modiJus  of  the  phase  1 

P2 

1.98  •  10® 

Pa 

sheEur  modulus  of  the  phase  2 

Ai 

94  •  10® 

Pa 

adiabatic  Lame  Lambda  coefficient  of  the  phase  1 

A2 

7.54 . 10® 

Pa 

adiabatic  Lame  Lambda  coefficient  of  the  phase  2 

cW 

1000 

J/(kg-K) 

specific  heat  capacity  of  the  phase  1 

c(^) 

2000 

J/(kg-K) 

specific  heat  capacity  of  the  phase  2 

aW 

2 • 10-® 

K-' 

volumetric  thermal  expansion  modulus  (phase  1) 

a(2) 

o 

1 

K-i 

volumetric  thermal  expansion  modulus  (phcise  2) 

Tl 

900 

s 

relaxation  time  of  the  phase  1 

T2 

10 

s 

relaixation  time  of  the  phcise  2 

fi 

25 

s-1 

Omega  (tensile  damage)  kinetic  factor 

C 

25 

s-' 

Alpha  (shear  deimage)  kinetic  factor 

D 

25 

s“^ 

Omega-Delta  (delamination  damage)  kinetic  factor 

A 

3.10® 

Pa*s 

Omega  enthropy  factor 

A 

3- 10® 

Pa*s 

Alpha  enthropy  factor 

Aa 

3-10® 

Pa*s 

Omega-Delta  enthropy  factor 

£. 

;  0.01 

- 

criticed  strain  factor 

0.01 

— 

criticed  shear  factor 

A. 

0.005 

- 

criticed  phases  fitness  factor 

/ 

10®  ^  10® 

Pa/s 

twisting  process  stress  per  time  factor 

/ 

10®  ^  10^° 

Pa/s 

straining  process  stress  per  time  factor 

Table  5.1.  Parameters  for  the  model  composite 


The  results  for  the  rate  /  =  10®  Pa/s  are  illustrated  in  Fig.  5.1(a-d).  The  strain 
eoz  grows  slowly  (Fig.  5.1a)  until  damage  parameter  wa  starts  accumulating  (Fig.  5.1c).  - 

Then  e0z  starts  growing  much  faster.  The  corresponding  time  <*,  strain  and  stress 

lb 


=  ft*  can  be  easily  detected  since  from  t  =  t*  the  behavior  of  £«*(<)  curve  differs 
sharply  from  that  in  the  absence  of  delamination  (t*;^  =  0,  see  the  dashed  curve  in 
Fig.  5.1a).  The  other  damage  parameter  a  remains  zero  under  these  loading  conditions 
(Fig.  5.1b).  The  temperature  variations  AT  are  shown  in  Fig.  5.1d.  While  all  the  d2miage 
pcirameters  are  equal  to  zero  AT  turns  to  be  slightly  negative  in  deforming  but  those 
negative  values  cannot  be  distinguished  within  the  scale  of  the  figure.  But  on  accumulation 
of  the  damage  parameters  in  the  present  case)  the  temperature  grows  rapidly  in 
irreversible  transformations  (Fig.  5.1d). 

The  Figs.  5.2a-d  illustrate  the  behavior  of  paraimeters  under  the  rate  of  loading  /  = 
10®  Pa/s.  It  is  seen  that  in  faster  loading  the  growth  of  starts  earlier  but  brings  to 
similar  picture  (Fig.  5.2a).  The  dashed  curve  on  the  Fig.  5.2a  shows  the  behavior  of 
egz  in  the  absence  of  delamination.  Damages  in  shear  have  not  occured  yet  (Fig.  5.2b). 
Temperature  vari^lrions  (Fig.  5.2d)  are  higher  tmder  these  loading  conditions. 

Under  even  faster  loading  (Fig.  5.3;  /  =  10^  Pa/s)  the  picture  is  qualitatively  the 
same.  Only  the  accumulation  of  damages  in  delcimination  starts  much  earlier  and  the 
criticsJ  vEdues  of  and  are  higher  for  the  present  case  thEin  that  for  the  previous 
cases.  This  result  shows  that  under  higher  rates  of  loading  the  critical  breaikup  limits  turn 
to  be  higher. 

Thus  we  can  have  a  number  of  experiments  under  low  rates  of  loading  to  determine  the 
critical  values  for  developing  damage  parcimeters  and  validation  of  results.  The 

results  for  temperature  variations  can  be  also  used  for  developing  the  model  parameters 
Eind  validating  purposes. 

The  increase  of  the  loading  rate  (/  =  10*  Pa/s)  radically  changes  the  scenarium 
of  the  process  (Fig.  5.4).  The  rapid  growth  of  strain  e$z  starts  much  earlier  due  to 
accumulation  of  damages  in  sheEur  (Fig.  5.4a);  the  dashed  curve  illustrates  the  strain  esz 
behavior  in  the  absence  of  sheEir  damages  accumulation.  The  dEunage  pEirameter  a  grows 
very  rapidly  (Fig.  5.4b)  while  the  damage  parEuneter  wa  remEiins  zero  (Fig.  5.4c).  The 
damages  in  delamination  have  not  enough  time  to  proceed  since  damages  in  sheEur  bring 
to  breakup  very  fastly.  The  temperature  growth  (Fig.  5.4d)  tEikes  place  due  to  irreversible 
transformations  in  sheEtr  damaging. 

The  further  increase  of  loading  rate  (/  =  10*  Pa/ s)  does  not  change  the  scenarium 
qualitatively  (Fig.  5.5).  The  criticEd  values  ej*  smd  =  /  •  t  for  the  present  cEise  csm  be 
determined  eis  well  Eind  they  turn  to  be  higher  thEin  that  for  the  previous  CEise  (Fig.  5.4a). 
Those  criticEtl  values  enable  to  develop  the  damage  parameters  responsible  for  damages 
in  sheEir  solving  the  inverse  problem  of  developing  model  parameters  bEised  on  the-results 
of  experiments. 
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Fig.  5.2  Twisting  of  a  tubular  sample  at  f  =  10®  Pa  /  s . 
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Fig.  5.3  Twisting  of  a  tubular  sample  at  f  =  10^  Pa  /  s . 


Fig.  5.4  Twisting  of  a  tubular  sample  at  f  =  10*  Pa  /  s 


5.2  Tension  of  a  tubular  sample 

In  pure  tension  of  samples  the  only  component  of  the  stress  tensor  different  from  zero  is 
assumed  to  change  lineeurly: 

=  f 't.  (5.2) 

The  experimentally  measured  values  <r„(t)  and  T{t)  depend  on  all  the  damage  parameters 
under  these  conditions:  damages  in  tension  a;,  in  shear  a  2ind  delamination  wa. 

The  low  rate  of  loading  {/  =  10®  Pa/s)  brings  to  the  case  of  slow  long-lasting  growth 
of  £**  (Fig.  5.6a)  that  is  changed  for  a  rapid  growth  when  delamination  occurs  (Fig.  5.6c). 
The  dashed  curve  on  the  Fig.  5.6a  shows  the  behavior  of  e**  in  the  absence  of  delamination. 
The  other  damage  parameters  remEiin  constant  (Figs.  5.6b, d).  Temperature  variation  in 
irreversible  transformations  is  relatively  not  very  high  (Fig.  5.6e).  The  critical  values  of 
e*  a*  (t*)  =  f  •  t*  can  be  determined  from  the  picture.  In  solving  the  inverse  problems 
of  developing  the  damage  constants  from  experiments  the  values  <r*,,  can  be  used 
for  validating  the  damage  in  delamination  constants  determined  in  the  experiments  on 
twisting  the  sample. 

The  increase  of  the  rate  of  loading  (/  =  10®  Pa/s)  does  not  change  the  picture  qualita¬ 
tively,  but'  quantitatively  decreases  the  breakup  time  and  increases  the  critical  'values 
and  (Fig.  5.7).  The  damage  par£imeters  a  and  a;  remain  equal  to  zero  (Figs.  5.7b,d) 
since  the  critical  values  and  are  still  much  lower  than  that  necessary  for  damages 
in  tension  and  shear  to  proceed. 

The  increase  of  loading  rate  up  to  /  =  10’^  Pa/s  does  not  change  the  scenarium 
qualitatively  (Fig.  5.8)  but  brings  to  essential  quantitative  changes.  The  critical  values 

and  <r*2  grow  up  and  approach  ones  under  which  the  damages  in  shear  and  tension 
can  occur. 

The  future  increase  of  the  loading  rate  (/  =  10®  Pa/s)  brings  to  the  scenarium  when 
breakup  is  achieved  due  to  growth  of  damages  in  shear  (Fig.  5.9b)  and  tension  (Fig.  5.9d). 
The  rapid  growth  of  e**  followed  by  the  breakup  (Fig.  5.9a)  takes  place  very  quickly  so 
that  dcimage  in  delamination  cannot  occur  (its  critical  parameters  and  <r*^  under  the 
present  loading  conditions  jure  higher  than  that  for  damages  in  shear  (e**  and  <r**)  and 
tension  {e***  and  aH*). 

The  dashed  curve  on  the  Fig.  5.9a  marked  ”a  =  0”  shows  the  behavior  of  £**  in  the 
absence  of  sheeur  damages  accumulation;  it  is  placed  very  near  to  the  actual  Szz  curve 
(solid  line),  that  means  that  the  breakup  under  present  loading  conditions  occurs  mostly 
because  of  tensile  damages  accumulation,  and  the  accumulation  of  damages  in  shear 
practically  does  not  influence  the  breakup  scenarium.  The  dashed  curve  on  the  Fig.  5.9a 
marked  ”u>  =  0”  shows  the  behavior  of  e**  in  the  absence  of  tensile  damages;  it  separates 
from  the  dashed  curve  marked  ”0!  =  a;  =  0”  later.  Comparison  of  the  dashed  curves 
marked  ”0;  =  0”  on  the  Figs.  5.9b,c,d  shows  that  in  case  of  neglecting  the  tensile  damages 
accumulation  the  breakup  will  be  delayed,  and  both  shear  and  delamination  will  take 
place.  However  the  delamination  does  not  take  place  when  we  account  for  the  the  tensile 
damages. 

Damages  in  tension  occur  a  little  bit  earlier  than  in  shear  (Fig.  5.9d)  and  their  growth 
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brings  to  strong  changes  in  e„(t)  function  behavior.  At  the  same  time,  the  successive 
growth  of  the  a-parameter  does  not  bring  essential  changes  to  the  growth  of  the  w- 
parameter  (Fig.  5.9d),  Ezz  strain  and  temperature  under  the  present  loading  conditions. 
The  dashed  curve  marked  "a  =  0”  in  Fig.  5.9d  shows  the  growth  of  u>  under  the  condition 
of  OL  (shear  damages  accumulation)  remaining  zero. 

The  critical  strains  e***  and  stress  c*”  can  be  determined  from  Fig.  5.9a. 

In  solving  the  inverse  problem  the  critical  values  determined  from  the  experimental 
curve  eiz(t)  could  be  used  for  developing  the  model  constants  for  damages  in  tension  (fl, 
A,  £*)  since  the  model  parameters  for  damages  in  shear  could  be  determined  from  the 
experiments  on  twisting  the  samples  (section  5.1). 

Fig.  5.10  illustrates  the  results  for  the  higher  rate  of  loading  /  =  10^  Pa/s.  It  is  seen 
that  accumulations  of  both  damage  parameters  a  and  u  contribute  to  the  rapid  growth 
of  Ezz  after  the  critical  values  are  reached  (Fig.  5.10a).  The  absence  of  one  of  the  damage 
parameters  a  or  a;  essentially  influences  the  growth  o^'^her  (Figs.  5.10b, d).  The  growth 
of  temperature  under  present  rate  of  loading  takes  place  mostly  due  to  the  growth  O'f  ot 
(damages  in  sheair). 

The  further  increase  of  the  rate  of  loading  (/  =  10^®  Pa/s)  does  not  change  the  picture 
quEilitatively  but  brings  to  an  increaise  of  the  critical  values  of  e***  and  a***  =  f-t***  and  to 
a  decrease  of  characteristic  breakup  time.  The  results  for  the  rate  of  loading  /  =  10^®  Pa/s 
Eire  shown  on  the  Fig.  5.11a-e. 
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Fig.  5.8  Tension  of  a  tubular  sample  at  f  =  10^  Pa  /  s . 


5.3 


Recommendations  on  developing  the  model  dam¬ 
age  constants  in  experiments 

The  numerical  simulations  of  dynamical  behavior  of  parameters  in  twisting  and  tension 
of  tubular  composite  samples  shown  definite  peculiarities  of  scenario  for  the  materials 
posessing  the  properties  close  to  that  given  in  the  table  5.1. 

For  the  low  rates  of  loading  both  in  twisting  and  tension  the  breakup  takes  place 
due  to  accumulation  of  damages  in  delamination  'while  the  damage  parameters  in  shear 
and  tension  remain  equal  to  zero.  Thus  the  critical  values  for  parameters 
ran  be  determined  from  the  experimental  strain-time  diagram  enabling  to  develop  the 
damage  constants  D,  Aa,  A*.  The  experimental  results  under  different  rates  of  loading 
ran  be  used  for  validating  purposes. 

For  the  high  rates  of  loading  in  twisting  the  bre2ikup  takes  place  due  to  accumuation 
of  damages  in  shear,  and  for  the  loading  in  tension  -  due  to  accumulation  of  damages 
both  in  she^u:  and  tension,  while  the  deimage  parameter  wa  remains  equal  to  zero.  Thus 
the  critical  values  ej*,  <rg*  can  be  determined  from  the  experiments  in  twisting.  Different 
rates  of  loading  can  provide  with  different  pairs  of  ej*,  (tJ*  that  would  make  it  possible 
to  determine  the  damage  constants  for  the  damage  in  shear:  C,  A,  £*.  The  experimental 
data  on  temperature  growth  AT  can  provide  with  additional  data  for  developing  those 
const^lnts  smd  validation. 

Experiments  in  high  rate  tension  can  provide  with  the  critical  parameters 
cheiracterizing  the  beginning  of  damaging  in  tension.  Assuming  the  model  constamts  for 
the  sheaur  damages  known  from  independent  experiments  in  twisting  one  can  develop  the 
model  constants  fl,  A,e*  for  the  damages  in  tension  baised  on  the  experimental  £**(<)  and 
T{t)  curves. 
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Chapter  6 

On  the  methods  of  modeling  the 
material  behavior  after  destruction 
criteria  having  been  satisfied  in 
some  zones 


Intensive  dyneimic  loading  of  the  structure  elements  (in  impact,  explosion  or  thermal 
loading)  brings  to  destruction  of  materials.  Dynamical  breaikup  is  a  complex  multistage 
process  involving  origination,  growth  and  coalescence  of  microdefects,  formation  of  mi¬ 
crocracks,  their  growth  to  macroscale,  intersections  leading  to  fragmentation  of  material. 

Two  characteristic  stages  can  be  distinguished  within  the  process  of  dynamical  breakup. 
The  first  stage  -  predestruction  (continual,  or  uniformly  distributed  destruction)  :  orig¬ 
ination  of  microdefects,  their  growth  and  coalescence,  formation  of  initial  macrocracks 
Eind  pores  in  the  process  of  irreversible  deformation  of  material. 

The  second  stage  -  destruction  of  material  :  formation  of  macrocracks  and  growth  of 
pores  inside  the  material,  their  coming  to  free  surface,  sepEiration  of  fragments. 

Mathematical  modeling  of  the  predestruction  stage  can  be  performed  introducing  dam¬ 
age  parameters  for  materials.  Models  of  this  type  applied  for  two-phase  laminated  com¬ 
posite  materials  were  described  in  [1,2]  and  are  discussed  in  the  chapters  4  cind  5  of  the 
present  report.  After  the  destruction  criteria  for  the  critical  specific  dissipation  is  satisfied 
in  some  place 

D  =  D.,  (6.1) 

the  present  computation  cell  gives  birth  to  a  new  destruction  interface  free  of  loadings 
(in  case  the  destruction  criterion  (6.1)  was  not  satisfied  in  the  adjacent  cells  before). 
Otherwise  the  crack  born  in  the  adjacent  cell  shotdd  grow  through  the  present  cell. 

Here  we  regard  the  approaches  to  mathematical  modeling  of  the  second  stage  of  de¬ 
struction.  There  exist  at  least  two  different  approaches  [3-6].  The  first  one  is  beised  on 
the  explicit  tracking  of  the  free  surfaces  of  the  crack.  The  second  is  replacing  the  de- 
structed  materiad  by  a  number  of  discrete  particles  characterized  by  their  sizes,  mass  eind 
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momentum.  Those  particles  can  interact  with  each  other  and  with  the  free  surfaces  of 

the  continuous  material  following  a  given  law. 

The  second  approach  is  not  in  fact  an  alternative  to  the  first  one  but  rather  its  further 
development.  The  second  approach  permits  to  create  a  simplified  procedure  for  numerics: 
the  destructed  material  is  assumed  to  be  a  media  resisting  compression  but  irresistible  for 
expansion  in  tension.  In  a  number  of  cases  this  simplification  is  significant  for  practical 
applications. 

Further  we  regard  the  both  approaches  taking  a  two-dimensional  case  as  an  example. 


6.1  The  method  of  explicit  free  surface  tracking  in 
macroscopical  destruction  of  material 

Using  the  Lagrangian  approach  for  description  of  the  media  deformation  the  cells  move 
with  the  media.  Thus  the  most  effective  method  of  free  surface  tracking  is  based  on  the 
local  reconstruction  of  the  grid  in  the  zone  of  the  crack  [3-6].  On  giving  birth  to  a  crack 
the  cell  is  removed  from  the  grid  and  replaced  by  the  two  free  surfaces  on  the  both  sides 
of  the  gap.  The  mass,  momentum  and  other  characteristics  of  the  cell  are  redistributed 
among  the  adjacent  ones.  The  both  sides  of  the  crack  are  either  free  of  tensions  or  have 
the  "contact”  boundeiry  conditions  depending  on  the  process.  If  the  crack  is  expanding 
a  free  surface  boundary  conditions  are  applied.  In  case  of  the  collapse  of  the  gap  the 
calculations  jure  conducted  using  the  contact  boundary  conditions  cdgorithm  [7]. 

The  basic  ideas  of  the  method  will  be  discussed  using  an  example  of  modeling  the 
border  sTirfaces  of  an  jurbitrary  curved  crack  in  a  continuum. 

Let  the  conditions  for  macro  destruction  (6.1)  be  satisfied  in  some  cell  of  the  grid.  The 
first  alternative  to  be  checked  under  conditions  is  either  birth  of  a  new  crack  takes  place 
or  the  further  propagation  of  a  crack  existing  in  the  adjacent  cell.  In  other  words,  we 
check  if  jiU  the  adjacent  cells  are  not  destructed. 

I.  In  case  of  a  new  crack  is  being  born  in  the  cell  the  following  algorithm  is  applied: 

1)  The  center  if  the  cell  where  the  destruction  criterion  (6.1)  was  satisfied  is  deter¬ 
mined: 


fc  =  1,2, 3,4  -  numbers  of  the  grid  points  surrounding  the  cell  (Fig.  6.1).  The  procedure 
is  performed  in  the  laboratory  system  of  coordinates  {®i,a;2}- 

2)  The  orientation  of  the  gap  plane  in  the  center  of  the  cell  is  determined.  The  orien¬ 
tations  of  the  planes  of  maximal  orthogonal  Cn  and  tangential  (7^  stresses  are  determined 
for  the  purpose  as  well  as  their  values  within  the  cell: 

O'n  =  (Tiia^ +  2(Ti20[fi  +  <r22fi^] 

=  ^(o-u«  +  0-12^)^  +  Wl20l  +  <^22^)*  - 


a*+^*  =  l; 

where  »  =  (a,/3)  -  a  unit  normal  vector  to  the  plane  of  maximal  stresses  (Fig.  6.1). 
The  procedure  described  above  enables  to  determine  two  planes  characterized  by  normal 
vectors  =  (a„,)3„)  and  =  {oLr.fir)^  The  first  one  (n„)  is  that  of  maximal  normal 
stress  (Tn  in  the  cell,  the  second  one  (ur)  is  the  plane  of  maximal  tangential  stresses  So 
one  has  an  jdternative  of  two  probable  orientations  of  the  originating  crack.  The  following 
criterion  permits  to  make  the  choice  of  the  direction  of  the  crack: 


max 


where  tb  is  the  meiximal  shear  limit,  as  -  meiximal  tensile  limit  (material  constants). 
If  the  condition 


is  satisfied  it  is  assumed  that  the  shear  crack  is  bom  within  the  cell  along  the  plane 
characterized  by  the  normsd  vector 
If  the  condition 


is  satisfied  within  the  cell  the  tensile  crack  is  to  be  born  along  the  plane  characterized  by 
the  normal  vector  n„.  The  criterion  formulated  above  is  the  classical  one  of  Davidenkov- 
Friedmann. 

3)  The  intersections  of  the  crack  plane  with  the  grid  hnes  surrounding  the  cell  are 
determined  (Fig.  6.1).  Let  the  intersection  A  be  on  (jfi,j2)  aJid  the  intersection  B  be  on 
(i3,i4)- 

4)  The  grid  points  and  are  moved  into  the  point  A,  the  grid  points  jsjji  are 
moved  into  the  point  B  (Fig.  6.2).  All  the  parameters  related  to  the  grid  points  and  the 
adjacent  cells  are  relocated  following  the  procedures  described  in  [8].  The  idea  of  the  grid 
modification  of  this  type  is  illustrated  in  the  Fig.  6.2. 


II.  In  case  of  propagation  of  the  crack  from  an  adjacent  cell  a  diflTerent  procedure  is 
used.  The  Fig.  6.2  shows  that  the  crack  can  propagate  further  either  in  the  cell  ii  or  in 
the  cell  12-  There  can  take  place  more  comphcated  cases  that  will  be  regarded  later. 

A.  Let  the  crack  grow  into  the  cell  ii,  i.e.  the  destruction  criterion  for  the  critical 
vedue  of  dissipation  (6.1)  is  satisfied  in  the  cell  *i.  And  the  values  of  the  normal  vector 
to  the  plane  n  =  are  determined.  Then  the  border  interfaces  of  the  crack  in  the  ti 

cell  can  be  determined  in  the  following  way. 

1)  Knowing  the  orientation  of  the  crack  plane  in  the  cell  ii  one  can  construct  the  plane 
beginning  from  point  C  of  its  intersection  with  the  border  of  the  cell  js,  je  shown  in 
the  Fig.  6.3. 

2)  The  grid  points  and  je  are  moved  to  the  point  C.  All  the  values  of  parameters, 
mass,  momentum,  energy  in  the  cells  is  and  are  redetermined  following  the  new  grid 
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pattern.  The  Fig.  6.4  presents  a  fragment  of  the  computational  grid  with  a  gap  in  two 
cells. 

B.  Let  the  crack  due  to  determined  orientation  of  maximal  stresses  planes  grow  outside 
the  cell  *1  (Fig.  6.2)  and  cross  the  cell  ia  for  example  (Fig.  6.3).  Then  tracking  the  borders 
of  the  crack  is  performed  by  the  following  algorithm. 

1)  The  gap  plane  is  constructed  beginning  in  the  point  A  and  tracked  until  it  intersects 
with  the  borders  of  the  cell  1*3.  Two  cases  are  possible: 

a)  the  point  of  intersection  D  is  found  on  the  border  J5J7  or 

b)  the  point  of  intersection  E  is  found  on  the  border  jVjg  (Fig-  6.5). 

2a)  In  case  ’’a”  the  grid  points  js,  jr  are  moved  into  the  point  D,  the  grid  points  Jq 
and  js  are  moved  into  the  center  of  AD.  The  new  grid  fragment  after  the  reconstruction 
is  shown  in  Fig.  6.6. 

2b)  In  the  opposite  case  ”b”  the  grid  points  jV  and  js  are  moved  into  the  point  E,  the 
grid  point  J5  is  moved  into  the  center  of  j^A  as  shown  in  Fig.  6.7. 

Thus  successive  satisfying  the  destruction  criterion  in  the  cells  enables  to  track  uniquely 
the  macroscopicjd  crack.  In  case  the  macroscopical  crack  comes  to  a  free  surface,  separa¬ 
tion  of  some  part  of  material  as  a  fragment  is  possible. 

The  precision  of  the  described  above  procedure  of  localization  of  a  crack  is  one  half  a 
mesh  width. 


6.2  Replacing  the  destructed  material  with  descrete 
particles 


Most  of  numerical  algorithms  are  constructed  in  a  way  providing  at  each  timestep  the 
following  data: 

1)  in  every  cell:  density  of  material,  stress,  strain  eind  strain  rates  tensors,  specific 
internal  energy,  dissipation,  temperature,  damadge  parameters  (in  case  the  model  for 
damadgeable  media  is  used); 

2)  in  the  grid  points:  velocities  and  current  coordinates.  Thus  the  first  group  of 
p^lrameters  is  related  to  the  centers  of  the  cells,  the  second  -  to  the  grid  points. 

In  case  the  destruction  criterion  (6.1)  is  satisfied  in  a  cell  the  material  in  the  cell  is 
considered  to  be  macroscopicaJly  damadged.  To  model  the  behavior  of  the  material  within 
the  damadged  cell  it  is  cissumed  that  the  material  resists  the  compression  only  and  has 
no  resistance  to  expansion  and  shecir.  That  mectns  the  stress  tensor  for  the  damadged 
medium  dij  is  a  spherical  one:  aij  =  —pSij.  The  equation  of  state  for  pressure  p  within 
the  damadged  ceU  has  the  following  form,  for  example: 


P=<i 


P>  P* 
P<P* 


(6.2) 


where  p*  is  the  density  of  the  materieil  within  the  ceU  by  the  time  of  destruction  (the  time 
for  the  criterion  (6.1)  is  satisfied),  K,n  -  material  constants. 
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In  case  the  damadged  ceU  is  an  internal  one,  the  calculations  are  run  following  the 
same  algorithm  but  for  the  equation  of  state  that  needs  to  be  taken  in  the  form  (6.2). 

In  case  the  damadged  cell  is  at  the  boundary  of  the  computational  domain  the  material 
of  the  cell  is  replaced  by  a  number  of  discrete  particles  with  the  radii  determined  froin  the 
condition  fitting  the  cell.  The  mass  of  the  cell  is  redistributed  among  the  particles.  Only 
one  layer  of  the  boundary  ceUs  can  be  transformed  into  discrete  particles  during  one  time 
step  because  it  is  assumed  that  the  velocity  of  the  destruction  wave  does  not  exceed  the 

velocity  of  acoustic  disturbances  in  the  media.  ,  .  „  ,  i,  -4.1, 

Thus  applying  such  an  algorithm  can  lead  in  some  cases  to  replacing  all  the  cells  witn 

discrete  particles,  i.e.  converting  a  continuous  body  into  a  cloud  of  fragments. 

Velocity  vectors  for  the  mass  points  (grid  points  and  discrete  particles)  are  determine 
from  the  following  finite-differences  equations; 

T  2tpi  ’ 


where  i  is  the  number  of  a  grid  point  or  a  particle,  n  -  the  number  of  a  time  step,  - 
mass  of  the  respective  mass  point,  r  -  time  step  value,  -  vector  of  forces  caused  by 
internal  stresses,  Ri  -  vector  of  reaction  forces,  that  is  equal  to  zero  in  the  mternal  grid 
points  and  is  determined  or  assigned  for  the  boundary  grid  points.  .111. 

Let  us  regard  the  most  general  case  when  the  computational  domain  contains  both  the 
grid  points  and  discrete  particles.  The  discrete  particles  are  assumed  to  be  incompressible 
and  interacting  with  the  boundaries  of  the  continues  medium  and  with  each  other.  An 
elementary  act  of  interaction  of  a  fragment  with  the  boundary  of  continuous  media  is 


modeled  in  the  following  way.  .  ^ 

The  known  by  the  time  <  =  t”  values  of  coordinate  vectors  xj*  and  velocities 
({  =  a,  6,  c)  for  the  grid  points  "a”,  "b”  and  the  discrete  fragment  "c”  make  it  possible  to 
determine  intermediate  values  of  coordinates  and  velocities  Xf  ,  (FlS-  6.8)  not  talung 
into  account  reaction  forces.  If  the  following  condition  is  satisfied  for  the  intermediate 

values  f,  e 

^Jabc  ^ 

<  Tc, 


h  = 


lab 


then  the  penetration  of  the  particle  "c”  into  the  continuous  material  must  l^ave  J;aken 
place.  Then  the  corection  of  coordinate  vectors  and  velocities  for  the  grid  points  ”a  ,  b 
and  the  particle  ”c”  is  necessary.  Here  is  the  particle’s  radius,  lab  -  the  distance  between 
”a”  and  "b”  grid  points,  h  -  the  height  of  the  triangle  "abc”.  Under  the  condition  h>rc 
there  is  no  interaction  of  particle  "c”  with  the  boundaries  of  the  continuous  material. 

To  perform  the  corrections  the  random  vector  components  being  the  result  of  mterac- 
tions  are  determined.  The  friction  forces  are  determined  as  well  if  the  law  is  given.  The 
grid  points  ”a”,  ”b”  and  the  particle  ’’c”  are  affected  by  the  following  forces: 

Ra  =  Nan  +  TcT,  Ra  = -{I  -  ot)Rc,  ^  (6-3) 
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C=|s.-ii,l,  a^rjr^,  (s,ri  =  o. 

Here  the  point  with  coordinates  xj,  is  the  projection  of  the  center  of  masses 

of  the  pmticle  onto  the  boundeiry  ah. 

The  value  of  the  tangential  force  Tc  is  determined  in  (6.3)  under  the  assumption  of  a 
no-slip  condition.  In  case  a  law  for  the  friction  is  given  the  value  Tc  determined  in  (6.3) 
remains  valid  when  the  following  inequality  is  satisfied: 

|Te|  <  k\Nc\,  (6.4) 

In  case  (6.4)  is  not  satisfied,  Tc  is  determined  by  the  formtila 

Tc  =  fc|i\r,|sign(r;),  (6.5) 

where  the  value  T*  is  given  by  (6.3).  The  value  of  k  in  (6.4),  (6.5)  is  a  friction  coefficient. 
It  is  evident  that  on  assuming  k  =  0  one  obtains  a  free  slip  condition  on  the  boundmy. 

Since  the  sizes  of  the  particles  are  final  (not  infinitely  small)  the  particles  interact  not 
only  with  boundaries  but  with  one  another  as  well.  Each  pmticle  is  chmacterized  by  its 
radius,  mass,  velocity  and  coordinate  of  its  center  of  masses. 

Let  the  line  segment  Uj  connect  the  centers  of  i-th  and  j-th  particles  (radii  and  Xj). 
If  the  condition 

lij  ^  Ti  r  j 

is  satisfied  then  the  particles  are  in  contact.  To  determine  the  reaction  forces  it  is  necessmy 
to  adopt  some  law  for  interaction  of  particles.  Using  the  analogy  with  the  interaction 
between  particles  and  the  boundciry  one  can  choose  the  law  of  absolutely  unelcistic  impact. 
Then  one  has  the  following  sets  of  formulas: 


S  __  ^li  -  % 
-  — - > 


rm, 


Riji  =  -Rlij ,  rn, 


1  1 

- 1 - 5 

rrii  rrij 


Uj  — 


X4 


(6.6) 


where  -  intermediate  value  for  the  vector  of  the  z-th  particle  center;  vi-  -  the  projection 
of  the  velocity  vector  on  the  fine  connecting  the  centers  of  peurticles;  Uj  -  a  unit  vector 
along  that  line. 

If  there  is  no  friction  between  the  pcirticles  then  the  tangential  component  of  the 
reaction  force  vector  is  zero,  and  there  is  no  necessity  in  the  corrections  of  the  tangential 
component  of  the  velocity  vector. 

If  there  exists  friction  given  by  the  law  (6.4)  then  computations  me  run  similar  to  the 
case  of  particle’s  interaction  with  a  boundmy  of  a  rigid  body  described  above. 

If  the  z-th  particle  interacts  with  a  number  of  pmticles  then  ciU  the  vectors  of  forces 
are  summed  up. 

If  one  adopts  any  other  law  of  pmticles  interactions  (that  is  possible)  then  it  is  neces¬ 
sary  to  define  an  appropriate  function  for  the  kinetic  energy. 
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6.3  Determining  the  number  of  fragments  within  a 
destructed  cell 

The  number  of  fragments  in  destruction  of  a  cell  after  the  dissipation  in  the  cell  reaches 
its  critical  value  cannot  be  arbitrary.  The  crack  formation  is  a  process  that  needs  energy. 
Thus  the  number  of  fragments  obtained  in  breakup  should  be  found  accounting  for  the 
balance  of  the  elastic  energy  accumulated  by  the  material  by  the  time  of  breakup,  and 
the  energy  necessary  to  form  all  the  free  surfaces  (cracks)  bounding  fragments: 


(6.7) 


where  S  is  the  total  area  of  the  free  surfaces  originating  within  the  cell,  7  -  energy 
necessary  for  breach  square  unit  formation,  E  -  elastic  energy  accumulated  within  the 
cell,  fee  -  coefficient  of  elastic  energy  transformation.  The  equation  (6.7)  still  leaves  some 
degrees  of  freedom  for  determining  the  number  and  masses  of  particles.  Thus  empirical 
data  on  particles  sizes  distributions  in  breakups  should  be  taken  into  account. 
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